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Abstract 

Intrinsic Delaunay triangulation (IDT) is a fundamental data structure in 
computational geometry and computer graphics. However, except for some 
theoretical results, such as existence and uniqueness, little progress has been 
made towards computing IDT on simplicial surfaces. To date the only way 
for constructing IDTs is the edge-flipping algorithm, which iteratively flips the 
non-Delaunay edge to be locally Delaunay. Although the algorithm is con¬ 
ceptually simple and guarantees to stop in finite steps, it has no known time 
complexity. Moreover, the edge-flipping algorithm may produce non-regular 
triangulations, which contain self-loops and/or faces with only two edges. In 
this paper, we propose a new method for constructing IDT on manifold trian¬ 
gle meshes. Based on the duality of geodesic Voronoi diagrams, our method 
can guarantee the resultant IDTs are regular. Our method has a theoretical 
worst-case time complexity 0(n 2 logn) for a mesh with n vertices. We ob¬ 
serve that most real-world models are far from their Delaunay triangulations, 
thus, the edge-flipping algorithm takes many iterations to fix the non-Delaunay 
edges. In contrast, our method is non-iterative and insensitive to the number 
of non-Delaunay edges. Empirically, it runs in linear time O(n) on real-world 
models. 

As a by-product, the regular Delaunay triangulations naturally induce dis¬ 
crete Laplace-Beltrami operators (LBOs), which are intrinsic to the geometry 
and have non-negative weights. We evaluate the commonly used discrete LBOs 
on the original triangulations and the intrinsic Delaunay triangulations. Com¬ 
putational results show that the IDT induced LBOs are more accurate than 
the LBOs defined on the original mesh. Moreover, their discrete Laplacian 
matrices have smaller condition number than the original triangulations. As a 
result, IDTs are ideal for applications which solve the linear system and eigen- 
system of the discrete Laplacian. 

Keywords: Intrinsic Delaunay triangulation, regular triangulation, geodesic 
Voronoi diagram, duality, discrete Laplace-Beltrami operator 
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1 Introduction 


A Delaunay triangulation for a set P of points in M 2 is a triangulation such that no point 
in P is inside the circumcircle of any triangle in the triangulation. It is well known that 
Delaunay triangulations tend to avoid skinny triangles, since they maximize the mini¬ 
mum angle of all the angles of the triangles in the triangulation. Although the Delaunay 
triangulations in Euclidean spaces are well understood [1], the intrinsic Delaunay trian¬ 
gulations on Riemannian manifolds have received less attention. By using the closed ball 
property [2], Dyer et al. |3] proposed adaptive sampling criteria for constructing intrinsic 
Voronoi diagram and its dual Delaunay triangulation on 2-manifolds. Recently, Boisson- 
nat et al. [5| proposed an algorithm for constructing intrinsic Delaunay triangulation on 
smooth closed submanifold of Euclidean space. Both methods are based on the convex 
neighborhood , which, in general, is an extremely small region around a point on the sur¬ 
face. As a result, in spite of their important theoretical values, they are not practical for 
piecewise linear surfaces, which are dominant in digital geometry processing. 

Rivin |5] and Indermitte et al. [6] defined intrinsic Delaunay triangulation (IDT) on 
triangle meshes, where the IDT edges are geodesic paths and the geodesic circumcircles 
of all Delaunay triangles have empty interiors. The intrinsic Delaunay triangulation has 
many nice properties. For example, Bobenko and Springborn [7] proved that the classic 
cotangent Laplace-Beltrami operator (LBO) has non-negative weights Wij , if and only if 
the underlying triangulation is Delaunay. They also proposed a new LBO which depends 
only on the intrinsic geometry of the surface and its edge weights are non-negative. 

Using intrinsic properties, such as the strong convexity radius and the injectivity radius, 
Dyer et al. [3] presented adaptive sampling criteria which guarantee the IDT is regular. 
Their elegant results establish inequalities that relate these intrinsic properties to the local 
feature size. However, to our knowledge, there is no practical algorithm for computing 
strongly convex regions on meshes. To date, the only practical algorithm for computing 
IDT is the edge-flipping algorithm [71IEJ 6], which iteratively flips the non-locally Delaunay 
edge to be locally Delaunay. The edge-flipping algorithm is conceptually simple and easy 
to implement. Indermitte et al. [6] showed that the edge flipping algorithm terminates 
in a finite number of steps and thus the intrinsic Delaunay tessellation exists. Bobenko 
and Springborn |7j proved the uniqueness of the intrinsic Delaunay tessellation. Although 
the edge flipping algorithm converges, it has no known time complexity. Moreover, the 
resultant IDT may contain non-regular triangles, which have either self-loops (i.e., edges 
with identical end points) or faces with only two edges. 
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In this paper, we propose a new method for constructing intrinsic Delaunay triangu¬ 
lation on 2-manifold triangle meshes. Our idea is to compute a geodesic Voronoi diagram 
(GVD), which has a regular dual triangulation. Given a 2-manifold mesh with n vertices, 
our method first computes the GVD by taking all vertices as sites. Then, our method 
identifies the Voronoi cells that violate the closed ball property [2] and fixes them by 
adding auxiliary sites. We show that by adding at most 0(n ) sites, the dual graph of 
the geodesic Voronoi diagram is an intrinsic Delaunay triangulation, which is guaranteed 
to be regular. Moreover, thanks to the bounded time complexity of computing geodesic 
Voronoi diagrams [9], our method has a theoretical worst-case time complexity 0{n 2 logn) 

. We observe that most real-world models are far from their Delaunay triangulations, thus, 
the edge-flipping algorithm takes many iterations to converge, whereas our method is not 
sensitive to the number of non-Delaunay edges and it empirically runs in linear time 0(n) 
on these models. 

As a by-product, the regular Delaunay triangulations naturally induce discrete Laplace- 
Beltrami operators (LBOs), which are intrinsic to the geometry and have non-negative 
weights. We evaluate the commonly used discrete LBOs on the original triangulations and 
the intrinsic Delaunay triangulations. Computational results show that the IDT induced 
LBOs are more accurate than the LBOs defined on the original mesh. Moreover, their 
discrete Laplacian matrices have smaller condition number than the original triangulations. 
As a result, IDTs are ideal for applications which solve the linear system and eigensystem 
of the discrete Laplacian. We demonstrate the IDTs on harmonic mapping and manifold 
harmonics, and observe that the IDTs produce more accurate and robust results than the 
original triangulations. 

The remaining of the paper is organized as follows: Section [2] reviews the mathematical 
background and highlights the fundamental difference between the 2D Delaunay triangula¬ 
tion and the IDT on meshes. Section [3] details our algorithm for computing IDTs, followed 
by the experimental results and comparison in Section [4j Section [5] shows the IDT induced 
discrete LBOs are more accurate and robust than the original meshes, and demonstrates 
them on solving the linear system and eigensystem of the discrete Laplacian. Finally, 
Section [6] concludes the paper. To ease reading, we list the main notations in Table [l] and 
delay the long proofs in Appendix. 
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M 

a manifold triangle mesh 




V, E, F 

the vertex, edge and face sets of M 




v,Vi,Vj, • • • 

vertex 




6, G-ij = Vj ‘ ‘ 

edge 




II 

-se 

triangular face 




p,q,Pi,qi,--- 

point on M 




lip, q) 

the geodesic path between p and q 




dip, q) 

the geodesic distance between p and q 



n{= \V\) 

number of vertices 




GVD(P) 

the geodesic Voronoi diagram of sites P 


T 

the set of Voronoi vertices 




B 

the set of Voronoi edges 




C 

the set of Voronoi cells 




VCipi) 

the Voronoi cell of site pi 




Kp, q) 

the bisector of sites p and q 




pKp) 

the pseudo bisector of site p in Section 

3.3 


V'(D v ) 

the set of augmented sites in Section 

3.3 


V"(D V') 

the set of augmented sites in Section 

3.4 


V"'(2 V") 

the set of augmented sites in Section 

3.6 


IDT{M ) 

the intrinsic Delaunay triangulation of M 



the set of g-edges 




r 

the set of geodesic triangles 




£ 

(/-edge 




T 

geodesic triangle 





Table 1: Main notations. 
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Figure 1 : Geodesic paths & geodesic triangles, (a) D\ is a simply connected domain with 
non-positive Gaussian curvature everywhere. There is a unique geodesic path between any 
pair of distinct points in D\. (b) D 2 is also simply connected. However, it has a part with 
positive Gaussian curvature. As a result, there are two geodesic paths between p and q. 
(c) D% has negative Gaussian curvature everywhere. However, as its fundamental group 
is non-trivial, the geodesic paths are not unique, (d) The red region r = (pi,P 2 ,P 3 ) is a 
geodesic triangle. However, its complement M \ r, although having three geodesic sides, 
is not a geodesic triangle, since it is of genus 1. (e) y(p 4 , £> 5 ), 7(ps, Pe) and 7(p6,P4) are 
three geodesic paths wrapping the cylinder. Each colored region has three geodesic sides. 
However, none of them is a geodesic triangle, since it is multiply connected. 

2 Preliminaries 

Let M be a manifold triangle mesh, and V, E , F be the set of vertices, edges and faces 
of M. Every interior point of M has a neighborhood which is isometric to either a neigh¬ 
borhood of the Euclidean plane or a neighborhood of the apex of a Euclidean cone. 

2.1 Geodesic Paths, Geodesic Triangles and Geodesic Circumcircles 

Consider two points p, q G M. A geodesic path between p and g, denoted by y(p, q), 
is the locally shortest path between them. Mitchell et al. m showed that the general 
form of a geodesic path is an alternating sequence of vertices and (possibly empty) edge 
sequences such that the unfolded image of the path along any edge sequence is a straight 
line segment and the angle of 7 passing through a vertex is greater than or equal to tt. 
The vertices with cone angles more than 2n are called saddle vertices , which play a critical 
role in geodesic computation nu. We denote by d(p, q) the geodesic distance between p 
and g, and 6 (p, q) the bisector of p and q. 

Given a simply connected domain Q with negative Gaussian curvature everywhere, 
the geodesic path y(p, q) is unique for any pair of points p, q G O. However, in general, 
geodesics are not unique for regions with positive Gaussian curvature and/or non-trivial 
topology. See Figure [lja)-(c) . 


Definition 1 (Geodesic Triangle). A geodesic triangle r G M is a simply connected do¬ 
main whose boundary dr has three geodesic paths. Each geodesic path is called a g-edge 
and the endpoints of a g-edge are called g-vertices. 


5 









Back view 


Front view Back view Front view 


Back view Front view 


Back view 


Front view 



Pi 
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(a) Apip 2 P 3 has no circumcircle (b) Ap^pQ has 2 circumcircles 


Figure 2 : Geodesic circumcircles. Geodesic triangle Apip 2 p^ has no geodesic circumcircle, 
since the bisector b(p\,p 2 ) (shown in blue) does not meet b(p 2 ,P 3 )- I n contrast, geodesic 
triangle Ap^pQ has two geodesic circumcircles, since b(p 4 , £> 5 ) (blue) and b(p§,p§) (red) 
intersect twice. 


On R 2 , any three non-parallel lines form a triangle. However, not all three geodesics 
on a mesh form a triangle. See Figure [ljd)-(e). 

Definition 2 (Geodesic Disk & Geodesic Circumcircle). A geodesic disk centered at a 
point p G M of radius r, denoted by D(p,r), consists of all points whose distance to p does 
not exceed r, i.e., D(p,r ) = {q G M : d(p,q) < r}. If a geodesic disk is simply connected, 
its boundary dD(p,r ) is called geodesic circumcircle. 

On R 2 , each non-degenerate triangle has a unique circumcircle that passes through its 
three vertices. However, a geodesic triangle may have no geodesic circumcircle at all or 
more than 1 geodesic circumcircles. See Figure [ 2 ] 

2.2 Intrinsic Delaunay Triangulations 

In contrast to Euclidean spaces, Delaunay triangulations do not exist for an arbitrary 
set of points on a Riemannian manifold. Leibon and Letscher m proposed sampling 
density conditions to ensure that the triangulation can accurately represent both the 
topology and geometry of the manifold. However, Boissonnat et al. [4] pointed out that 
sampling density alone is insufficient to guarantee an intrinsic Delaunay triangulation for 
Riemannian manifolds of dimension 3 and higher. 

The definition of intrinsic Delaunay triangulation on 2-manifold meshes is due to 
Rivin [5], who generalized the 2D Delaunay condition (i.e., a circle circumscribing any 
Delaunay triangle does not contain any other input points in its interior) and required the 
empty geodesic circumcircle property. 
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Definition 3 (Intrinsic Delaunay Triangulation). The intrinsic Delaunay triangulation 
(IDT) on M, denoted by IDT(M) = (V, 5,T), is a triangulation such that 

• the vertex set of IDT(M) equals V; 

• every edge £ in S is a geodesic path on M, i.e., a g-edge; 

• each face r G T is a geodesic triangle, which has a geodesic circumcircle containing 
no mesh vertices in its interior; 

• and T forms a tessellation of M; 


Indermitte et al. [6] proved the existence by showing the edge-flipping algorithm ter¬ 
minates in finite steps. Bobenko and Springborn [Tj proved the uniqueness of Delaunay 
tessellation (whose faces are generally but not always triangular). The Delaunay triangu¬ 
lation can be obtained by triangulating the non-triangular faces. They pointed out that 
the Delaunay triangulation, while in general not unique, differs from another Delaunay 
triangulation only by edges with vanishing cot-weights (i.e., the sum of two angles opposite 
that edge is 7r). 

2.3 Geodesic Voronoi Diagrams 

Definition 4 (Geodesic Voronoi Diagram). Let P = {pi,P2? • • • ,Pm} be a set of points 
on M. The Voronoi cell VC(pi) corresponding to site pi consists of all points whose dis¬ 
tance to pi is less than or equal to their distance to any other site, i.e., VC(pi) = {q G 
M\d(pi,q) < d(pj,q),\/i ^ j}. The geodesic Voronoi diagram (GVD) of P is the union of 
all Voronoi cells, GVD(P) = {VC(pi), VC(p 2 ), • • • , VC(p m )}. The Voronoi edges bound¬ 
ing the Voronoi cells are trimmed bisectors and the Voronoi vertices are points incident to 
three or more Voronoi edges. 

It is shown in mm that the GVD forms a tessellation of M, since all Voronoi cells 
are mutually exclusive, and \Jqf 1 VC(pi ) = M. 

It is well known that Voronoi cells in M 2 are simply connected and convex, and the 
Voronoi edges are all line segments. However, these properties do not hold for geodesic 
Voronoi diagrams. Although a geodesic Voronoi cell is still connected, it may have multiple 
boundaries and even handles. See Figure [3] (right). Moreover, a geodesic Voronoi cell 
can be concave, since the bisectors on triangle meshes consist of both line segments and 
hyperbolic segments |9j. See Figure [3] (left). 
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Figure 3: Bisectors and Voronoi cells on a 2-manifold mesh are significantly different from 
their Euclidean counterpart. Left: There is a saddle vertex (yellow) on the geodesic path 
y(p, q) (cyan). As a result, the bisector b(p, q ) (in blue) consists of both line segments and 
hyperbolic segments. Right: VC(ps ) (in green) is of genus-1, and VC (pi) (in red) and 
VC(p2) (in yellow) are of genus-0 but with multiple boundaries. 


Voronoi diagram Bisectors Dual graph 



Figure 4: Consider the Two-hole Torus model with 8 sites p^ i = 1, • • • , 8. Note that 
VC(p§), VC(pq), VC(pr) and VC(ps) are multiply connected. VC(pi) and VC(p 2 ) share 
2 edges, and so do the other pairs of Voronoi cells (VC(p 2 ),VC(ps)), (VC (ps) ,VC (pa)) , 
(VC(pi), VC(p 4 )) and (VC (P 2 ) ,VC (p 4 )). We can clearly see that the dual graph is not a 
triangulation. In fact, it does not even form a tessellation. 


2.4 The Closed Ball Property & Duality 

It is well known that the dual graph of a Voronoi diagram in M 2 is the Delaunay tessellation 
for the same set of points. Therefore, one can adopt the Voronoi diagram’s algorithm to 
construct the Delaunay tessellation, and vice versa. This duality, however, does not exist 
on 2-manifold meshes in general. See Figure [4j 

Edelsbrunner and Shah [2j introduced the closed ball property for triangulating an 
abstract topological space. For a surface without boundary, the closed ball property 
expresses three conditions: 

1. Homeomorphism condition: each Voronoi cell is homeomorphic to a planar disk; 

2. 2-cell intersection condition: the intersection of any two Voronoi cells is either 
empty, or a single Voronoi edge, or a single Voronoi vertex; 

3. 3-cell intersection condition: the intersection of any three Voronoi cells is either 
empty or a single Voronoi vertex. 

Dyer et al. [3] gave the sufficient condition that a GVD has a dual triangulation. 
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Algorithm 1 Constructing Intrinsic Delaunay Triangulation from the Dual of Geodesic 
Voronoi Diagram 

Input: M — (V,E,F), a 2-manifold triangle mesh 
Output: IDT(M ), the intrinsic Delaunay triangulation of M 
1: GVD(V) — compute_gvd(M) (Procedure [2] in Section [372] ) 

2: GVD(V ') = ensure Jiomeomorphismmondition(GVD(V)), where V C V' (Procedure 
[3] in Section 3.3) 

3: GVD(V") = ensme42-cellJntersectionmondition(GVD(V')), where V' C V" (Proce¬ 
dure [4] in Section 3.4) 

4: IDT(M) = compute dual graph(GyD(y // )) (Procedure [ 5 ] in Section 3.5) 


Theorem 1 |3j. If a geodesic Voronoi diagram satisfies the closed ball property, its dual 
IDT exists. Moreover, the IDT is regular so that 1) each geodesic triangle has three distinct 
g-edges and three distinct g-vertices; and 2) any two geodesic triangles, if their intersection 
is not empty, share either a common g-vertex or a common g-edge. 


Dyer et al. also proved that if there are at least four distinct sites in the GVD and 
both the homeomorphism condition and 2-cell intersection conditions are satisfied, then 
the 3-cell intersection condition is redundant. 


3 Constructing Intrinsic Delaunay Triangulations 


This section presents the algorithm for constructing intrinsic Delaunay triangulation on 
meshes. In Sections 3.1-3.5, we assume the input 2-manifold mesh M is closed. We then 


consider open meshes in Section 3.6, In Section 3.7, we analyze the time complexity of 
our algorithm. 

We also assume that the meshes are free of degenerate cases, i.e., no four vertices lying 
on a geodesic circumcircle. These degenerate cases can be easily handled by the symbolic 
perturbation technique m ■ 

3.1 Overview 


Flipping edges and computing the dual of Voronoi diagrams are two commonly used tech¬ 
niques for constructing Delaunay triangulations in R 2 . As mentioned before, the edge-flip 
algorithm on meshes PS does not have a known time complexity and it may produce 
non-regular triangulations, which have self-loops or faces with only two edges. In this 
paper, we take the other direction by computing the dual of geodesic Voronoi diagrams. 
The major challenge in this direction is that not every GVD has a dual Delaunay triangu- 
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(a) GVD{V ), V = {vijlh (b) GVD(V '), V' = ^U(<7i}(c) GVD(V") : V" = V' UW (d) Dual of GVD(V") 
Figure 5: Algorithmic pipeline, (a) Using the original vertices as sites, the geodesic Voronoi 
diagram does not have a regular dual Delaunay triangulation, since the Voronoi cell VC(v i) 
is not simply connected, (b) Adding a new site q\ to V ensures the homeomorphism 
condition so that both VC{yi) and VC(qi) become topological disks. However, they 
share two common edges (see the front and back views), (c) Adding another site s, which 
is equidistant to v\ and #i, the new geodesic Voronoi diagram GVD(V") satisfies the 
closed ball property, (d) The dual graph of GVD{V") is a regular Delaunay triangulation. 
Color schemes: the mesh vertices are red dots and the added auxiliary sites are green. The 
mesh edges are grey, the Voronoi edges blue, the Delaunay edges green, and the geodesics 
purple, respectively. 


lation. As shown in Section |T4} the closed ball property is the sufficient condition for the 
existence of a dual Delaunay triangulation. Therefore, our goal is to enforce the closed 
ball property everywhere on the GVD. We consider only the homeomorphism condition 
and the 2-cell intersection condition in our algorithm, since all the models we are dealing 
with have more than 4 vertices. 

Our algorithm (ref. Algorithm [TJ consists of four steps. Firstly, taking all mesh 
vertices as sites, it computes the geodesic Voronoi diagram GVD(V). Secondly, it checks 
the homeomorphism condition for all Voronoi cells. If a Voronoi cell, say VC(vi), is 
not homeomorphic to a disk, the algorithm adds an auxiliary site in VC{yi) and then 
locally updates the GVD. Thirdly, it checks the 2-cell intersection condition for all pairs 
of adjacent Voronoi cells. If two adjacent Voronoi cells, say VC{vi) and VC(vj), have more 
than 1 common Voronoi edges, the algorithm adds an auxiliary site that is equidistant 
to Vi and Vj, and locally updates the GVD. After steps 2 and 3, the updated GVD is 
guaranteed to have a dual triangulation, which is computed in Step 4. 

Intuitively speaking, our method adaptively increases the sampling density for the 
regions where the closed ball property fails. Figure [5] illustrates the algorithmic pipeline 
using a toy cylinder model. 
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3.2 Computing the GVD 


The first step in our algorithm is to construct the geodesic Voronoi diagram on M. Liu 
et al. [9] presented a generic algorithm, which runs in 0{n 2 \ogn) time for m(< n) sites, 
where n = \V\ is the number of vertices in M. One, of course, can apply Liu et al.’s 
algorithm directly to our application. However, our scenario is slightly different in that 
we take all mesh vertices as sites. Moreover, as our goal is to compute the dual Delaunay 
triangulation, we don’t need to explicitly construct the Voronoi diagram. Therefore, we 
adapt Liu et al’s algorithm for IDT construction (see Procedure [2]). We will show in 


Section [377] that our adapted GVD algorithm runs empirically in linear time on real-world 
models. 

Procedure 2 compute_gvd(M) 

1: Taking {^i,^ 2 5 * * • ,v n } as sources, compute the geodesic distance by the MMP algo¬ 
rithm m- 

2: Identify the triangles at which three or more bisectors meet. 

3: Compute the symbolic representation of Voronoi edges and Voronoi cells. 
Postcondition: The symbolic representation of GVD(V ) 


Using the mesh vertices as sources, we compute the geodesic distances by the Mitchell- 
Mount-Papadimitriou (MMP) algorithm [15]. Each mesh vertex is assigned 0 (since it is 
the source) and each mesh edge is partitioned into disjoint intervals, called windows , where 
a window encodes the shortest distance to some source vertex. Edge e contains a bisector 
if e has windows corresponding to different sources. Then, we identify the triangles where 
three or more bisectors meet. Rather than computing the geometry of each Voronoi cell 
explicitly, we need only a symbolic representation: the bisector of and vj is an ordered 
pair {vi,Vj} with i < j, and a Voronoi cell is an ordered list VC{vi) — {v ^, V { 2 , • • • }, where 
each pair of vi and vi k corresponds to a trimmed bisector, i.e., a Voronoi edge. 

In case the geometry of a Voronoi cell is required, which is rarely encountered in IDT 
construction on real-world models, the bisectors can be computed on-the-fly. To compute 
the boundary edge b(vi,Vi k ) of VC(vi), we find the triangles containing the two branch 
points {vi,Vi k _nVi k } and {vi,Vi k ,Vi k+1 }, then we locally trace the bisector by unfolding 
the corresponding triangles onto R 2 . 

3.3 Ensuring the Homeomorphism Condition 

In R 2 , the term “bisector” refers to the set of points which are equidistant to two distinct 
sites. In contrast to the 2D counterpart, there are two types of bisectors on a mesh: one 
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is the conventional bisector of two distinct sites, and the other is a special bisector, called 
pseudo bisector , which corresponds to a single site. 

Definition 5 (Pseudo-bisector). The pseudo-bisector of a site v consists of points q G 
VC(v ) such that there are two or more geodesic paths between v and q of equal length. 

Pseudo bisectors usually occur in a cylinder-shaped region. See Figure ga). The 
following proposition states that a pseudo bisector, if exists, must be a line segment after 
unfolding. 

Proposition 1. Let GVD(P) be the geodesic Voronoi diagram of sites P, where V C P. 
If, for a site p G P, the Voronoi cell VC(p), has a pseudo bisector pb(p), then pb(p) is a 
line segment when making all faces containing pb(p) coplanar. 

Since pseudo bisectors are due to low sampling density, we add an auxiliary site on 
each pseudo-bisector to destroy it. Given a Voronoi cell VC(pi) which contains a pseudo 
bisector pb(pi ), we compute a point q G pb(pi) such that q minimizes the geodesic distance 
d(pi,x ), Vx G pb(pi). If the unfolded images of j(pi,q) and pb(pi) are perpendicular, q 
is the intersection point. Otherwise, q is one of pb(pif s two end points. We then locally 
update the GVD around the new site q. The following proposition guarantees that Voronoi 
cell VC(pi ) becomes a topological disk when the auxiliary site q is added. 

Proposition 2. Assume that Voronoi cell VC(pi) has a pseudo-bisector pb(pi) and the 
point q G pb(pi) minimizes the geodesic distance d(pi,x) for all x G pb(pi). Add the auxil¬ 
iary site q into the GVD. Then both VC(q) and VC(pi) are topological disks. 


Procedure 3 ensure_homeomorphism_condition(GVD(V)) 

1: V' = V 

2: for every v t G V do 

3: if VC(vi) has a pseudo-bisector pb{vi) then 

4: Compute q G pb(vi) that minimizes d(vi,x ), x G pb(vi) 

5: V' = V{J{q} 

6: Locally update GVD(V f ) 

7: end if 

8 : end for 

Postcondition: Each Voronoi cell in GVD(V') is homeomorphic to a disk 
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Proposition 3. GVD(V') has at most An + 4g — 4 Voronoi vertices, 6n + 6g — 6 Voronoi 
edges, and 2n Voronoi cells, where n — \V\ and g is the genus of M . 



Figure 6: Voronoi cells VC{vi) and VC(vj ) violate the 2-cell intersection condition, since 
they share two Voronoi edges b & and ty. To destroy one Voronoi edge, say &&, we add an 
auxiliary site sG&fc and update the GVD locally. The gray shaded areas are some Voronoi 
cells, which are irrelevant to the discussion. 

3.4 Ensuring the 2-cell Intersection Condition 

After adding auxiliary sites on the pseudo-bisectors, all Voronoi cells in GVD(V') are 
simply connected. Let B denote the set of Voronoi edges (i.e., trimmed bisectors). Recall 
that the symbolic representation of the bisector of and vj is an ordered pair {vi,Vj} 
satisfying i < j. We sort all the Voronoi edges {vi,Vj} in B in the ascending order by 
their first elements. If two ordered pairs have the same first element, the second element 
is used to break a tie. 

If an order pair {vi,Vj} appears more than once in the sorted list B , the two corre¬ 
sponding Voronoi cells VC{vi) and VC(vj) violate the 2-cell intersection condition. Let 
bk and bi be the two common edges shared by VC(vi) and VC(vj). See Figure |6^a). To 
destroy one edge, say, b we compute a point s G b^ which minimizes the geodesic distance 
d(pi,x ), Vx G bk. As bk bisects pi and pj , the path 7 (pj, s) is a minimizing geodesic. We 
then add s into V' and update the Voronoi cells locally. See Figure j6^b). The following 
proposition guarantees that the updated Voronoi cells satisfy both the homeomorphism 
condition and the 2-cell intersection condition. 

Proposition 4. Let VC{vi) and VC(vj) be adjacent Voronoi cells that share two Voronoi 
edges b & and The point s G minimizes d{vi,x) for x G Then the GVD with 
additional site s satisfies both the homeomorphic condition and the 2-cell intersection con¬ 
dition. 
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Procedure 4 ensure_2-cell_intersection_condition(GPi?(P / )) 

1 : Let V" = V 7 

2 : Sort all Voronoi edges in the ascending order 

3: for any two Voronoi edges bk = {pi,Pj} and bi = {pi,Pj} with the same ordered pair 

do 

4: Compute point sG&fc that minimizes d(pi,x), x G b^ 

5: V" = V"{J{s} 

6 : Locally update GVD(V") 

7: end for 

Postcondition: GVD(V") satisfies the closed ball property 


Proposition 5. GVD(V") 


has 0{n ) Voronoi vertices, Voronoi edges and Voronoi cells. 



Figure 7: Each Voronoi edge corresponds to a unique Delaunay edge. Consider two 
adjacent Voronoi cells VC(p ) and VC(q). Let D = VC(p){jVC(q). Assume there are 
two geodesic paths y(p, q) G D and y'(p, q) G D, which connect p and q. Let T (resp. J 7 ') 
denote the set of faces that y(p, q) (resp. 7 '(p, g)) passes through. Since y(p, g) ^ 7 7 (p, g), 
the two face sets are different T 7 ^ Tj meaning that there is at least one vertex, say v, 
inside the region bounded by y(p, q) and q). As a result, the Voronoi cell VC(v) is in 
between VC(p ) and VC(q ), which is a contradiction. Therefore, there is a unique geodesic 
path y(p, q) G D between p and q , which is taken as the Delaunay edge dual to b(p , q). 


3.5 Computing the Dual Graph 

Since the closed ball property holds everywhere for GVD(V // ), the dual Delaunay trian¬ 
gulation exists. By the 2-cell intersection condition, two adjacent Voronoi cells VC(p ) and 
VC(q) share only one common Voronoi edge, denoted by b(p,q). As Figure [T| shows, the 
Voronoi edge b(p , q) corresponds to a unique Delaunay edge, which is dual to b(p , q). 

Given a Voronoi edge b(p,q) in GVD(V "), we compute its dual Delaunay edge by 
using the geodesic distance field obtained in Procedure 2 (line 1). Note that the MMP 
algorithm splits each mesh edge into two oriented halfedges with opposite directions and 
each halfedge contains the windows, a discrete data structure that encodes the geodesic 
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paths to the closest source. Since both VC(p) and VC(q ) do not contain any mesh vertex 
other than the generating sites p and q, each side of b(p , q ) contains exactly one window, 
which encodes the geodesic paths to p and q, respectively. We denote by T v (resp. F q ) the 
set of faces containing the geodesic paths from p (resp. q) to any point on b(p,q). Then 
the unique Delaunay edge y(p, q) is computed by unfolding the faces U T q . 


Procedure 5 compute_dual_graph(GPD('\/ // )) 

1: for each Voronoi edge b(p,q) G B in GVD(V") do 

2: Using the windows stored at the two sides of b(p, q), compute the face sets T v and 

T q , respectively. 

3: Compute y(p, q) by unfolding the triangles in T v |J Fq- 

4: end for 

5: for any three Voronoi cells that meet at a Voronoi vertex do 
6: Create a geodesic triangle with the three corresponding g-edges. 

7: end for 

Postcondition: IDT(M) is a regular IDT. 


(a) 



(b) 


Figure 8: The closed ball property for manifolds with boundary, (a) VC(pi) PidM has two 
non-adjacent boundary Voronoi edges (in green), which violate the condition Al. Case 
(b) VC(pi ) fl VC(pj) fl dM has two boundary Voronoi vertices (blue dots), which violate 
the condition A2. The boundary dM is shown in red. 



3.6 Manifolds with Boundaries 

Now we consider meshes with boundaries. Let us denote by dM the boundary of M. The 
closed ball property for a manifold with boundary has two more conditions [2]: 

Al. The intersection of any single Voronoi cell and dM is empty, a single point or a 
single line segment homeomorphic to [0,1]. 

A2. The intersection of two distinct Voronoi cells and dM, VC(pi) f) VC(pj) f) dM, i ^ 
j, is either empty or a single point. 

These additional conditions complement to the three conditions of the closed ball 
property for closed manifolds. Figure [8] shows two cases where either condition Al or A2 
does not hold. These issues can be fixed by adding auxiliary sites to the boundary dM. 


Proposition 6. The Voronoi cell VC{v) in GVD(V") has two non-adjacent boundary 
Voronoi edges b & and bi (i.e., bbi G dM and b^ Dbi = 0^. Without loss of generality, say 
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v £ bk. The point s G ^ minimizes the geodesic distance d(v,x) for x G bThe GVD 
with the additional site s satisfies both boundary conditions Al and A 2. 


Procedure 6 ensure_boundary_condition(Gy,D(y // )) 

1: Let V" 1 = V" 

2: if a Voronoi cell VC(v ) has two non-adjacent boundary Voronoi edges then 
3: Set bk be the boundary Voronoi edge that does not contain v 

4: Compute point s G bk that minimizes d(v,x), xG^ 

5: V'" = V'"\J{s} 

6: Locally update GVD(V f ") 

7: end if 

Postcondition: GVD(V"') satisfies the boundary conditions Al and A2 


Proposition 7. GVD(V"') has 0(n) Voronoi vertices, Voronoi edges and Voronoi cells. 

3.7 Complexity Analysis 

Theorem 2. Given a closed 2-manifold mesh M — (V,E,F), our IDT construction 
algorithm (Algorithm 1) has a worst-case time complexity 0{n 2 \ogn), where n— \V\. 

Proof. Our algorithm consists of four steps, corresponding to Procedures 2, 3, 4 and 5. 

First, we show that Step 1, computing the geodesic Voronoi diagram GVD(V ), takes 
0{n 2 \ogn) time, which is the dominant term of all four steps. The GVD algorithm is 
built upon the MMP algorithm m, a classic discrete geodesic algorithm which represents 
the discrete geodesic wavefront using windows. Maintaining the windows in a priority 
queue, the MMP algorithm iteratively propagates windows across the faces and updates 
the geodesic distance when a window covers a vertex or part of an edge. It terminates 
when the priority queue is empty, i.e., all vertices and edges have been covered by some 
windows. Mitchell et al. showed that each edge has at most 0(n) windows, thus, the 
MMP algorithm has a worst-case time complexity 0{n 2 \ogn). In the GVD computation, 
0(n) windows on an edge e means e has 0{n) bisectors. Then lines 2 and 3 in Procedure 

2 have 0(n 2 ) time complexity. Therefore, Step 1 takes 0(n 2 log n) time. 

Then we show that Steps 2, 3 and 4 take 0(n 2 ) time. Given a Voronoi cell VC(pi) that 
contains a pseudo-bisector pb(pi ), finding the point q that minimizes d(pi,x), Vx G pb(pi ), 
takes 0(1) time. Locally updating VC(q) takes at most 0(n) time. By Proposition 3, 
there are at most n auxiliary sites added in Step 2. Thus, Step 2 takes 0{n 2 ) time. 
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(b) Normal case 

Figure 9: Complexity of bisectors, (a) The triangles on the side of the cylinder are global , 
since they span the whole model. Each of these global triangle contains 0(n) bisectors. As 
a result, computing GVD(V ) takes 0(n 2 log n) time, (b) In general, if the majority of mesh 
faces are local (i.e., each triangle covers only a local region on M), most edges intersect 
only 1 or 2 bisectors (see the close-up view). Therefore, the GVD can be computed in 
0{n) time. 


Similarly, by Proposition 5, there are 0(n ) auxiliary sites added in Step 3, and updating 
the GVD for each new site takes 0(n) time. Therefore, Step 3 takes 0(n 2 ) time. 

In Step 4, computing a Delaunay edge y(p, q) takes 0{n) time, since there are at most 
0(n ) faces in the set T p |J T q . As GVD(V") has 0(n) Voronoi edges, Step 4 also takes 
0(n 2 ) time. 

Putting it all together, our algorithm has a worst-case time complexity 0(n 2 logn). 

□ 


The theoretical worst-case time complexity 0{n 2 logn) is very pessimistic, since it happens 
only when each triangle contains 0{n ) bisectors. See Figure [9](a) for a model with the 
worst-case time complexity. We call a triangle with 0(n) bisectors global , since it indeed 
spans a global region on the model. We observe that on many real-world models, the 
majority of the triangles are not global, even though the mesh triangulation is poor (i.e., 
far from its Delaunay triangulation). Computational results show that on average, a mesh 
edge on a real-world model has only 0(1) bisectors. See Figure^b). As a result, all of the 
four steps in our algorithm run in 0(n) time. Computational results in Figure [To] confirm 
that our algorithm has an empirical linear time complexity. 
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♦Eight 

♦Kitten 

♦Lucy 

♦Buddha 


Figure 10: Our algorithm runs in linear time empirically. The horizontal axis shows the 
mesh complexity and the vertical axis is the execution time (in seconds). 


Model (|E|) 

(&meam &maxi &std) 

Edge-flipping 

Ours 

cond(A) 

Time 

Nfe 

Time 

Nas 


M 

IDT(M) 

CSG (357) 

(4.74 

,33.02,2.13) 

0.006 

48.4% 

0.004 

1 

1.3 

X 

10 5 

4.9 

X 

10 :i 

Fandisk (IK) 

(44.2, 

177.5, 9.49) 

0.015 

37.3% 

0.011 

0 

2.4 

X 

10 7 

7.7 

X 

10 4 

Eight (3K) 

(1.67, 

11.80, 0.88) 

0.031 

23.6% 

0.016 

0 

7.6 

X 

10 4 

3.9 

X 

10 4 

Sphere (6K) 

(1.67, 

59.60, 2.14) 

0.019 

0.08% 

0.006 

0 

5.5 

X 

10 5 

5.5 

X 

10 5 

Teapot (17K) 

(1.69, 

20.74, 1.66) 

0.171 

23.5% 

0.125 

0 

8.7 

X 

10 5 

3.4 

X 

10 5 

Decocube (24K) 

(1.56, 

13.43, 0.65) 

0.188 

17.6% 

0.140 

0 

5.3 

X 

10 5 

2.5 

X 

10 5 

Fertility (37K) 

(2.34, 

17.20, 1.11) 

0.609 

46.3% 

0.561 

0 

8.5 

X 

10 8 

2.0 

X 

10 8 

Crank (60K) 

(17.5 

,279.5,14.2) 

1.539 

47.6% 

1.350 

0 

1.9 

X 

io 10 

9.5 

X 

10 6 

Bunny (216K) 

(1.23, 

11.83, 0.19) 

0.780 

1.40% 

0.156 

0 

4.1 

X 

10 6 

2.6 

X 

10 6 

Armadillo (519K) 

(1.31, 

95.29, 2.39) 

2.324 

4.59% 

0.983 

0 

1.1 

X 

10 7 

3.2 

X 

10 6 

Lucy (789K) 

(1.45, 

34.24, 1.38) 

5.179 

9.44% 

3.960 

0 

5.8 

X 

10 7 

4.3 

X 

10 7 

Buddha (1,196K) 

(1.47, 

29.16, 5.22) 

8.502 

8.66% 

5.379 

0 

6.8 

X 

10 7 

3.6 

X 

10 7 


Table 2: Mesh complexity and performance statistics. |E|: the number of edges in M; 
Nf e (%): percentage of the edges that are flipped by the edge-flipping algorithm; N as : 
the number of auxiliary sites added by our algorithm; cond(A): the condition number of 
the discrete Laplacian matrix. The 3-tuple (cr mean , cr max , &std) shows the mean, max and 
standard deviation of the anisotropy measure a. The running time was measured on an 
Intel Core i7-2600 CPU (3.40 GHz). 
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Figure 11: Results. The original mesh edges and the Delaunay edges are colored in black 
and red, respectively. The figures are of high resolution, allowing close-up examination. 
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4 Experimental Results 


This section reports the experimental results and compare the performance of our method 
and the edge-flipping algorithm. 

We implement our algorithm in C++ and test it on 10 synthetic and real-world models 
with diverse geometric and topological features. As Table [2] shows, most meshes are far 
from their Delaunay triangulations. For example, the Fertility model has 46.3% non- 


Delaunay edges (see Figure 11), which takes the edge-flipping algorithm many iterations 
to fix them. Our algorithm, in contrast, computes the IDT in a non-iterative manner and 
its performance is not sensitive to the number of non-Delaunay edges. We observe that 
our method consistently outperforms the edge-flipping algorithm in terms of execution 
time. 

We also investigate the relation between mesh quality and performance. We measure 
the quality of a triangle t by its anisotropy a — where p is the half-perimeter, h is 

the length of its longest edge and S is the triangle area. It is easy to verify that a(t) > 1 
and the equality holds when t is equilateral. We measure the triangulation quality by 


the mean 


l5 maximum a„ 


and standard deviation a st d of a for all triangles of M. 


Usually, a larger a means the higher degree of anisotropy and the further away the mesh 
is from its IDT, thus, the more edges that are flipped by the edge-flipping algorithm, and 
the higher the speedup our algorithm provides. 

As mentioned above, the edge-flipping algorithm may produce non-regular IDT. In 
contrast, our method guarantees the regular IDT by adding auxiliary sites at the poorly- 
sampled region. See Figures [l3| and 14 Theoretically, our algorithm adds 0(n) auxiliary 
sites to ensure the existence of dual triangulation. In practice, we observe that only a 
very small number of auxiliary sites are required for real-world models. This is not a 
surprise, since the auxiliary sites are only added on the sparsely-sampled regions which 
are not homeomorphic to a disk (e.g., the cylinder-like geometry in Figures [5] and 14). As 
Table [2] shows, although most test models are far from their Delaunay triangulations, they 
have a fairly good sampling density. Consequently, the resultant GVD(V ) automatically 
satisfies the closed ball density, and our algorithm does not add any new site at all. For 
these models, both our method and the edge-flipping algorithm produce exactly the same 
results. 

Our method can be easily adapted to centroidal Voronoi tessellation (CVT), which is a 
special type of Voronoi diagram such that the site of each Voronoi cell is also its center of 
mass. The existing work of CVT focuses on the convergence [16] > isotropic meshing m , 
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Figure 12: Centroidal Voronoi tessellation. Row 1: The existing CVT algorithms (e.g. d) 
produce non-regular Voronoi cells when the generating sites are sparse. The green arrow 
indicates a multiply-connected Voronoi cell and the pink arrows show Voronoi cells sharing 
2 common edges. Row 2: As a post-process, our method fixes the issues by adding auxiliary 
sites to the poorly-sampled regions and locally optimizing the site locations. Row 3: The 
dual graph of the updated CVT is a valid triangulation. The original sites and new sites 
are colored in red and green, respectively. 


energy functional smoothness m, and intrinsic computation m • However, none of them 
addresses the issue of regularity. Indeed, all of them simply take it for granted that the site 
number is sufficiently large so that the CVT is regular. Unfortunately, when the generating 
sites are sparse , non-regular Voronoi cells do exist. Figure [T2] shows an example of the 
genus-1 Kitten model with only 50 sites. Due to the low sampling density, the Voronoi cells 
around the tail do not satisfy the closed ball property. As a post-process to the existing 
techniques, our method automatically identifies the non-regular cells, adds auxiliary sites 
and then locally optimizes the site’s location using the Lloyd method PH- The updated 
CVT is regular, leading to a valid dual triangulation. 

5 Discrete Laplace-Beltrami Operators on Intrinsic Delau¬ 
nay Triangulation 

The Laplace-Beltrami operator (LBO) is an important operator on Riemannian manifolds 
and it has many desired properties. For example, applying the LBO to the coordinate 
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Figure 13: Comparison with the edge-flipping algorithm. Consider a prism with 7 vertices. 
Among the four edges facing V 2 V 3 and v^vq are local Delaunay, whereas the other 
two edges V 2 V 3 and v^vq violate the local Delaunay condition. Thus, the edge-flipping 
algorithm iteratively flips the non-Delaunay edges until the local Delaunay property holds 
everywhere. However, the resulting Delaunay triangulation is not regular, due to the self¬ 
loop (in red) connecting v\ and itself. Our method adds an auxiliary site v% which is the 
middle point of the self-loop and refines the Voronoi diagram locally. As a result, the dual 
Delaunay triangulation is regular. The bottom row shows the 2D flattening of the region 
in which the Delaunay property does not hold. 



Figure 14: The edge-flipping algorithm produces a self-loop on the CSG model (see the 
left and middle images). Our method solves the problem by adding an auxiliary site (see 
the green dot in the right image). Note that the original mesh edges (in blue) are line 
segments, whereas the IDT edges (in red) are geodesic paths. 


functions gives the mean curvature; the eigenfunctions of the Laplacian form a natural 
basis for square integrable functions on the manifold analogous to Fourier harmonics for 
functions on a circle. 

5.1 Convergence & Accuracy 

There is considerable amount of work of defining LBO on discrete domains 1201 EH E21 ESI 
EUZU2S1ES1EZU2S]. Most of these methods are variants of the cotangent scheme EH], 
which is a form of the finite element method applied to the Laplace operator on a surface: 

A /(^) = 52 w v - f( v j)) > 

{vi,vj}eE 
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where f : V M is a piecewise linear function defined on M. The weight for edge 

eij = {vi,vj} is 



( 1 ) 


where aij and otji are the two angles facing edge e^. The function / is discrete harmonic, 
if A f(vi) = 0 at all interior vertices. In spite of its extreme popularity in digital geometry 
processing, the cotangent formula has two serious drawbacks: 


The edge weight w\j in Equation ([TJ could be negative, meaning that fi is not a con¬ 


vex combination of its neighbors. In parameterization, this non-convex combination 
issue produces flipped triangles in the parametric domain. 

• The formula is not intrinsic, that is, two surfaces that are isometric but with 
different triangulations may have different discrete LBOs. 

Bobenko and Springborn [7] proved that the cotangent weight is non-negative if the un¬ 
derlying triangulation is Delaunay. Moreover, since the intrinsic Delaunay tessellation is 
unique, the discrete LBO defined on IDT(M) is also intrinsic and independent of the 
triangulation of M. 


In this subsection, we thoroughly evaluate the accuracy and convergence of commonly 
used discrete LBOs on the original mesh M and the intrinsic Delaunay triangulation 
IDT(M). 

(1) Classic cotangent LBO [20]: 


A (1) /Oi) = - f( v j)) 


{vi,Vj}eE 


where W{j is defined in Equation ([!]). 

(2) 1-ring-area-weighted cotangent LBO [22]: 


A <2| /(».) = P) /(«.) 


where A is the sum of areas of ^’s 1-ring triangles. 

(3) Voronoi-area-weighted cotangent LBO [|2'3 j : 


A(3) /(^) = A Wf( Vi ) 


where A vorono i is the Yoronoi area at vertex [29j. 

(4) Mesh LBO [2Sj : 


-voronoi 
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where A(f) is the area of triangle /. The parameter h is a positive quantity, which 
intuitively corresponds to the size of the neighborhood considered at Vi. 

We test the above discrete LBOs on R 2 and § 2 where the analytical LBO is available: 


d 2 f d 2 f 
A R 2/(X,2/) = g^ + ^2 


and 


A s2 /(M) = - 


1 d 


sin6% I + 


1 d 2 f 
30 J ' sin 2 6 dcf> 2 ’ 


sin 0 86 

where 6 and </> are the spherical coordinates. 

Following [26], we consider three functions on R 2 , fi(x,y) = x 2 , f 2 (x,y) = e x and 
/ 3 (x, y) = e^ +2/ , and three functions on 8 2 , y, z) = x 2 , /s(x, y, z) = and fe(x, y, z) = 
e ^+y < We generate a sequence of triangle meshes with increasing resolution for each func¬ 
tion /(#, y). To obtain the planar meshes, we apply the greedy triangulation [30] to 
random samples in [—0.5, 0.5] x [—0.5, 0.5]. We adopt the marching cube algorithm [31] to 
construct the spherical meshes, where the resolution is related to the user-specified cube 
size, i.e., the smaller cube size, the higher mesh resolution we obtain. As Table [3] shows, 
the IDTs induced LBOs produce smaller normalized L 2 error than those of the original 
meshes. The classic cotangent formula A^ does not converge to the analytical LBO at 
all, whereas the other three discrete LBOs converge. Moreover, when the mesh has a 
sufficiently high resolution, A^ 4 ) has the least L 2 error. 


5.2 Mean Curvature Computation 

We evaluate the performance of various discrete LBOs on mean curvature computation. 
Applying the LBO to the coordinate functions of a point p G S', we obtain the mean 
curvature vector, i.e., 

A p = 2H(p)n(p) 

where H(p) and n (p) are the mean curvature and unit outward normal at p. Following 
1231, we consider smooth surfaces given by the following non-linear functions: 

Fi(x,y) = ^4 - (x - 0.5) 2 - (y - 0.5) 2 , 


F 2 (x, y ) = tanh(9y - 9x), 

1.25 + cos(5.4y) 
F3(X ’ V) = 6 + 6(3s-1)2 ’ 

F 4 (x,y) = e -fi((*-o.5) 2 +0/-o.5) 2 ) 
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Planar domain (unit 

square 

[-0.5,0. 

5] x [— 

0.5,0.5]) 





Normalized L 2 error 

Function 

\E\ 

AW 

A^ 

AW 

AO) 



M 

IDT 

M 

IDT 

M 

IDT 

M 

IDT 


261 

0.731 

0.666 

0.655 

0.656 

0.905 

0.321 

2.801 

2.399 

h(x,y) 

1,121 

0.803 

0.844 

0.604 

0.595 

0.391 

0.183 

0.587 

0.513 

4,641 

0.960 

0.942 

0.548 

0.520 

0.362 

0.153 

0.186 

0.171 


18,881 

0.961 

0.960 

0.469 

0.417 

0.255 

0.119 

0.041 

0.040 


261 

0.839 

0.768 

0.641 

0.575 

0.745 

0.379 

6.626 

6.353 

-P (nr* 

1,121 

0.956 

0.879 

0.640 

0.571 

0.452 

0.160 

1.248 

1.058 

J 2 {x,y) 

4,641 

0.986 

0.883 

0.639 

0.569 

0.420 

0.149 

0.386 

0.334 


18,881 

0.996 

0.874 

0.632 

0.543 

0.334 

0.135 

0.082 

0.074 


261 

0.843 

0.741 

0.675 

0.623 

0.425 

0.230 

2.838 

2.786 

hix,y) 

1,121 

0.951 

0.849 

0.647 

0.621 

0.373 

0.161 

1.305 

1.225 

4,641 

0.988 

0.884 

0.644 

0.613 

0.257 

0.131 

0.332 

0.313 


18,881 

0.997 

0.924 

0.637 

0.578 

0.200 

0.105 

0.076 

0.067 



Spherical domain 

(unit sphere) 






Normalized L 2 error 

Function 

\E\ 

AW 

A^ 

Aw 

A^ 



M 

IDT 

M 

IDT 

M 

IDT 

M 

IDT 


804 

0.892 

0.892 

0.654 

0.654 

0.831 

0.300 

2.984 

3.053 

h(x,y,z) 

3,468 

0.972 

0.972 

0.642 

0.642 

0.432 

0.190 

0.561 

0.545 

14,268 

0.992 

0.992 

0.642 

0.642 

0.434 

0.177 

0.154 

0.148 


57,684 

0.998 

0.998 

0.641 

0.641 

0.368 

0.171 

0.027 

0.027 


804 

0.671 

0.697 

0.648 

0.663 

0.817 

0.394 

6.930 

6.804 

h (x,y,z) 

3,468 

0.881 

0.791 

0.604 

0.600 

0.408 

0.153 

1.222 

1.140 

14,268 

0.893 

0.877 

0.543 

0.536 

0.341 

0.119 

0.319 

0.330 


57,684 

0.932 

0.892 

0.442 

0.445 

0.248 

0.093 

0.054 

0.054 


804 

0.652 

0.701 

0.722 

0.699 

0.396 

0.211 

3.508 

2.858 

fe(x,y,z ) 

3,468 

0.813 

0.829 

0.604 

0.603 

0.393 

0.171 

1.275 

1.182 

14,268 

0.859 

0.878 

0.540 

0.555 

0.316 

0.156 

0.280 

0.268 


57,684 

0.881 

0.891 

0.431 

0.464 

0.301 

0.145 

0.053 

0.053 


Table 3: Evaluation of various discrete LBOs on planar and spherical domains, which 
are triangulated in increasing resolutions. In each resolution, the domain with the least 
normalized L 2 error is shown in red. \E\ is the number of edges in the meshes. 
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Function 

\E\ 

Normalized L2 error 

AW 

A^ 

A^ 


M 

IDT 

M 

IDT 

M 

IDT 

M 

IDT 


539 

0.861 

0.805 

0.580 

0.491 

0.081 

0.039 

2.378 

2.023 


2,311 

0.868 

0.819 

0.580 

0.443 

0.078 

0.033 

2.362 

1.765 

Fi(x,y) 

9,423 

0.869 

0.783 

0.580 

0.354 

0.078 

0.029 

2.324 

1.455 


38,177 

0.870 

0.820 

0.580 

0.319 

0.078 

0.024 

2.235 

1.317 


153,095 

0.870 

0.818 

0.580 

0.253 

0.078 

0.019 

2.004 

0.953 


334 

0.691 

0.655 

0.522 

0.469 

0.253 

0.186 

0.661 

0.535 


1,666 

0.698 

0.643 

0.481 

0.351 

0.080 

0.048 

0.654 

0.473 

F2 (x,y) 

7,372 

0.699 

0.657 

0.470 

0.310 

0.045 

0.024 

0.648 

0.447 


30,655 

0.700 

0.652 

0.468 

0.255 

0.039 

0.019 

0.635 

0.370 


124,431 

0.700 

0.673 

0.467 

0.192 

0.038 

0.014 

0.628 

0.270 


547 

0.752 

0.721 

0.519 

0.458 

0.064 

0.056 

0.854 

0.752 


2,677 

0.758 

0.691 

0.510 

0.391 

0.049 

0.028 

0.848 

0.606 

F3 (x,y) 

11,469 

0.759 

0.751 

0.508 

0.332 

0.039 

0.025 

0.832 

0.572 


47,437 

0.760 

0.718 

0.507 

0.291 

0.041 

0.020 

0.793 

0.447 


192,825 

0.760 

0.699 

0.507 

0.203 

0.041 

0.014 

0.683 

0.328 


401 

0.825 

0.752 

0.693 

0.553 

0.473 

0.401 

0.712 

0.633 


1,687 

0.829 

0.780 

0.684 

0.509 

0.462 

0.333 

0.700 

0.501 

F±{x,y) 

7,117 

0.830 

0.765 

0.682 

0.436 

0.461 

0.276 

0.675 

0.470 


29,265 

0.830 

0.800 

0.682 

0.362 

0.460 

0.252 

0.626 

0.338 


118,673 

0.830 

0.789 

0.682 

0.309 

0.459 

0.195 

0.553 

0.234 


Table 4: Mean curvature computation. Among the four discrete LBOs, the Voronoi-area- 
weighted LBO defined on the IDT (in red) produces the most accurate results. \E\ is the 
number of edges in the meshes. 


Similar to the convergence test, we generate a sequence of triangle meshes with in¬ 
creasing resolutions for each smooth surface F{x,y). Then we evaluate the accuracy of 
the mean curvature computed by using the conventional LBOs and the IDT induced LBOs. 
As Table [4] shows, the IDT induced LBOs produce more accurate results than the LBOs 
defined on the original meshes. 

5.3 Applications 

As mentioned above, the IDT induced cotangent LBO is intrinsic to the geometry and 
non-negative for all edges. These features are highly desirable to many graphics appli¬ 
cations, such as denoising [22] . parameterization [32j 33J, quadrangulation [341 . manifold 
harmonics [35>], shape signature [36] . diffusion distance m and biharmonic distance [38] . 
just name a few. In this subsection, we demonstrate IDT on harmonic mapping and 
manifold harmonics, two typical applications based on the discrete LBO. The conformal 
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discrete Laplacian matrix l c m is 


{ Xfc w iki if i=j 

-Wij, if {vi, Vj} G E 
0 otherwise 

where wij is the cotangent weight for edge { 77 , Vj} (see Equation 0). Let A be a diagonal 
matrix, where An is the Voronoi area at vertex Then the discrete Laplacian matrix 
L [23] is given by 

L = A -1 L C . 


The harmonic mapping is to solve a linear system of L. Let 0 : V -4 1 be a scalar 
function defined on mesh vertices. Since rank( L) = n — 1 , we need to fix at least one 
vertex to get a unique solution. Without loss of generality, say (j){y i) = 0. The function </> is 
harmonic if A= 0 for all vertices v^i> 1. The discrete harmonic map (j) is realized by 
solving the linear system L <fi = b, where </> = (..., 0(^),.. .) T , b = (1,0,0, • • • ) T G R nXl , 
and matrix L coincides with the discrete Laplacian matrix L except for the first row, which 
is (1, 0, 0, • • •). We compute the condition number of L, which is a good measure of the 
numerical stability of the linear system. As Table [i] shows, the IDT induced matrix L has 
smaller condition number than the matrix produced by the original mesh. 

We use harmonic mapping to parameterize the Fandisk model, which is a genus-0 
model with 1 boundary (i.e., a topological disk). The boundary vertices are mapped 
to unit circle S 1 using arc-length parameterization (see Figure |T5](b)) . We evaluate the 
parameterization quality by measuring the angle distortion 

E cot aa 2 + cot /3b 2 + cot 7 c 2 , 

- 73wT\ - A U )> 


feF 


4A (/) 


where a, 6, c, a , /?, 7 are the side lengths and angles of triangle / G F, /' G M 2 is the 
parameterized 2D triangle of /, and A(-) is the triangle area. To visualize </> on the IDT 
using texture mapping, we tessellate each non-planar geodesic triangle to planar polygons 
and linearly interpolate the texture coordinates for points at which the geodesic edges and 
the mesh edges meet. As Figure [15] shows, the IDT based harmonic mapping is numerically 
more stable than the original mesh. 

Manifold harmonics [35] are a natural generalization of the Fourier transform to curved 
domains. Since the conformal discrete Laplacian matrix L c is real symmetric, its eigen¬ 
values are real and its eigenvectors are orthogonal, which define a function basis allowing 
for such a transform. We transform the geometry into frequency space by projecting 
the coordinate functions onto the eigenvectors corresponding to low frequency. Then we 
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(a) Triangulation (b) 2D domain (c) Texture mapping 


Figure 15: We parameterize the Fandisk model to 2D unit disk using harmonic map. The 
original mesh (row 1) has many skinny triangles, so there are many negative weights in 
the cotangent LBO. Our IDT (row 2) guarantees non-negative weights for all internal 
edges and enforces the convex combination for all internal vertices. As a result, the IDT 
based paramterization is more stable and has less artifact than the original mesh. The 
parameterizations of M and IDT(M) have angle distortion 1.70 and 1.55, respectively. 


transform the object back into geometry space using the inverse transform. Figure 16 


shows the geometry reconstructed from manifold harmonics with an increasing number of 
eigenvectors. We observe that the IDT based manifold harmonics have less artifacts than 
those based on the original meshes, especially when only a small number of eigenvectors 
are used for reconstruction. 
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